The data being used in this workflow is taken from “Spatial Gene Expression Changes in the Mouse Heart After Base-Targeted Irradiation”
Radiation cardiotoxicity (RC) is a clinically significant adverse effect of treatment for patients with thoracic malignancies. Clinical studies in lung cancer have indicated that heart substructures are not uniformly radiosensitive, and that dose to the heart base drives RC.
This was a pilot experiment to determine whether spatial transcriptomics could provide useful insights in this context.
The data is included in the posit cloud workshop but if you are running this workshop elsewhere you can down load it from dropbox
A subset of the data is used for this workshop and it includes:
2 samples 2 conditions - irradiated heart and control 2 adjacent sections taken from each tissue block (sample) 4 sections in total
The sample metadata is summarised below:
dataDir <- "testData/results"
samples <- list.files(dataDir)
samples_short <- samples %>%
gsub("results_", "", .) %>%
gsub("_Results", "", .)
sample_id <- samples_short %>%
gsub("_.*", "", .)
condition <- c(rep("Control", 2), rep("irradiated", 2))
tissue <- c(c(rep("tissue1", 2), rep("tissue2", 2)))
sampleInfo <- data.frame(Fname = samples,
id = samples_short,
condition,
tissue,
sample_id
)
sampleInfo
## Fname id
## 1 results_C01_QUB_002_B1NP_Control_Results C01_QUB_002_B1NP_Control
## 2 results_C02_QUB_002_B1NP_Control_Results C02_QUB_002_B1NP_Control
## 3 results_C05_QUB_003_B2BP_irradiated_Results C05_QUB_003_B2BP_irradiated
## 4 results_C06_QUB_003_B2BP_irradiated_Results C06_QUB_003_B2BP_irradiated
## condition tissue sample_id
## 1 Control tissue1 C01
## 2 Control tissue1 C02
## 3 irradiated tissue2 C05
## 4 irradiated tissue2 C06
The data is exported from the 10X pipeline Spaceranger and it conists of the aligned and quantified reads for the transcriptomics data and the processed image data.
pathToTenXOuts <- file.path(dataDir, sampleInfo$Fname, "outs")
The SpatialFeatureExperiment library contains an import function for 10X data which takes the path to 10X files and loads the data in as a SpatialFeatureExperiment object
sfe <- SpatialFeatureExperiment::read10xVisiumSFE(
dirs = pathToTenXOuts,
samples = sampleInfo$Fname,
sample_id = sampleInfo$sample_id,
type = "sparse",
data = "filtered",
images = "hires",
zero.policy=TRUE)
saveRDS(sfe, "rObjects/sfe.rds")
The import function is a bit time consuming so save time I have pre-computed a loaded data rds file which can be imported here:
sfe <- readRDS("rObjects/sfe.rds")
sfe
## class: SpatialFeatureExperiment
## dim: 31053 14053
## metadata(0):
## assays(1): counts
## rownames(31053): ENSMUSG00000051951 ENSMUSG00000089699 ...
## ENSMUSG00000096730 ENSMUSG00000095742
## rowData names(1): symbol
## colnames(14053): AAACAAGTATCTCCCA-1 AAACAGAGCGACTCCT-1 ...
## TTGTTTGTATTACACG-1-3 TTGTTTGTGTAAATTC-1-1
## colData names(4): in_tissue array_row array_col sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
## imgData names(4): sample_id image_id data scaleFactor
##
## unit: full_res_image_pixel
## Geometries:
## colGeometries: spotPoly (POLYGON)
##
## Graphs:
## C01: col: visium
## C02: col: visium
## C05: col: visium
## C06: col: visium
Built of SpatialExperiment Class (which is built on SingleCellExperiment)
Incorporates geometries and geometry operations with the sf package
## SP and SF
sp and sf are to classes used in Spatial Analysis. sp is the older of the 2 and uses S4 objects, it has been around for 20 years but there is a move towards sf (simple feature) that sf uses S3 classes.
We currently use both sp and sf in our package but we are currently
working towards migrating to 100% sf workflow
https://github.com/r-spatial/sf/wiki/migrating
This function adds the spot positions to the as the spot coordinates which are used for plotting. For more details the functions are included R/sfe_prepareSFE.R
sfe <- Spaniel::prepareSFE(sfe, sampleInfo)
## Joining with `by = join_by(sample_id)`
Different regions of the heart were annotated using the imaging software Amira. These could also be created using Photoshop or other imaging software. Each annotation was then exported as a separate layer. The annotation images are the same size as the high resolution images exported from
For example this image shows the left ventricle in sample “C01”:
C01
Left Ventricle sample C01
The annotation_images directory
annotationDir <- "testData/annotation_images/"
annotations <- list.files(annotationDir)
annoInfo <- data.frame(sample_id = gsub("_.*", "", annotations) ,
domain = gsub(".*_", "", annotations) %>%
gsub("\\.jpg", "", .),jpg = annotations)
annoInfo
## sample_id domain jpg
## 1 C01 centre C01_centre.jpg
## 2 C01 left-atrium C01_left-atrium.jpg
## 3 C01 left-ventricle C01_left-ventricle.jpg
## 4 C01 right-atrium C01_right-atrium.jpg
## 5 C01 right-ventricle C01_right-ventricle.jpg
## 6 C01 valves C01_valves.jpg
## 7 C01 XR-LV C01_XR-LV.jpg
## 8 C01 XR-RV C01_XR-RV.jpg
## 9 C02 centre C02_centre.jpg
## 10 C02 left-atrium C02_left-atrium.jpg
## 11 C02 left-ventricle C02_left-ventricle.jpg
## 12 C02 right-atrium C02_right-atrium.jpg
## 13 C02 right-ventricle C02_right-ventricle.jpg
## 14 C02 valves C02_valves.jpg
## 15 C02 XR-LV C02_XR-LV.jpg
## 16 C02 XR-RV C02_XR-RV.jpg
## 17 C05 centre C05_centre.jpg
## 18 C05 left-atrium C05_left-atrium.jpg
## 19 C05 left-ventricle C05_left-ventricle.jpg
## 20 C05 right-ventricle C05_right-ventricle.jpg
## 21 C05 valves C05_valves.jpg
## 22 C05 XR-LV C05_XR-LV.jpg
## 23 C05 XR-RV C05_XR-RV.jpg
## 24 C06 centre C06_centre.jpg
## 25 C06 left-atrium C06_left-atrium.jpg
## 26 C06 left-ventricle C06_left-ventricle.jpg
## 27 C06 right-atrium C06_right-atrium.jpg
## 28 C06 right-ventricle C06_right-ventricle.jpg
## 29 C06 valves C06_valves.jpg
## 30 C06 XR-LV C06_XR-LV.jpg
## 31 C06 XR-RV C06_XR-RV.jpg
Radiation field = black rectangle; right atrium = yellow; left atrium = blue; great vessels = red; right ventricle = purple; left ventricle = green; atrioventricular valves = orange; central fibrous body = black; irradiated right ventricle = white; irradiated left ventricle = pink.
annoInfo$radiation_field <- ifelse(annoInfo$domain %in% c("left-ventricle",
"right-ventricle"),
yes = FALSE, no = TRUE)
annoInfo
## sample_id domain jpg radiation_field
## 1 C01 centre C01_centre.jpg TRUE
## 2 C01 left-atrium C01_left-atrium.jpg TRUE
## 3 C01 left-ventricle C01_left-ventricle.jpg FALSE
## 4 C01 right-atrium C01_right-atrium.jpg TRUE
## 5 C01 right-ventricle C01_right-ventricle.jpg FALSE
## 6 C01 valves C01_valves.jpg TRUE
## 7 C01 XR-LV C01_XR-LV.jpg TRUE
## 8 C01 XR-RV C01_XR-RV.jpg TRUE
## 9 C02 centre C02_centre.jpg TRUE
## 10 C02 left-atrium C02_left-atrium.jpg TRUE
## 11 C02 left-ventricle C02_left-ventricle.jpg FALSE
## 12 C02 right-atrium C02_right-atrium.jpg TRUE
## 13 C02 right-ventricle C02_right-ventricle.jpg FALSE
## 14 C02 valves C02_valves.jpg TRUE
## 15 C02 XR-LV C02_XR-LV.jpg TRUE
## 16 C02 XR-RV C02_XR-RV.jpg TRUE
## 17 C05 centre C05_centre.jpg TRUE
## 18 C05 left-atrium C05_left-atrium.jpg TRUE
## 19 C05 left-ventricle C05_left-ventricle.jpg FALSE
## 20 C05 right-ventricle C05_right-ventricle.jpg FALSE
## 21 C05 valves C05_valves.jpg TRUE
## 22 C05 XR-LV C05_XR-LV.jpg TRUE
## 23 C05 XR-RV C05_XR-RV.jpg TRUE
## 24 C06 centre C06_centre.jpg TRUE
## 25 C06 left-atrium C06_left-atrium.jpg TRUE
## 26 C06 left-ventricle C06_left-ventricle.jpg FALSE
## 27 C06 right-atrium C06_right-atrium.jpg TRUE
## 28 C06 right-ventricle C06_right-ventricle.jpg FALSE
## 29 C06 valves C06_valves.jpg TRUE
## 30 C06 XR-LV C06_XR-LV.jpg TRUE
## 31 C06 XR-RV C06_XR-RV.jpg TRUE
The annotations can be loaded in and added to the sfe object using the spaniel function findAllDomains. This creates a separate for each histological domain in the colData of the object. The combineDomains function can then be used to combine the annotations into a single column named domain.
These functions are found in: R/sfe_findDomain.R
For each of the images the functions:
sfe <- Spaniel::findAllDomains(sfe, annoInfo, annotationDir)
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## Warning in imager::grayscale(im): Image appears to already be in grayscale mode
## create a column containing all domains
sfe <- combineDomains(sfe, unique(annoInfo$domain))
sfe$domain
## add radiation field
sfe$radiation_field <- colData(sfe) %>%
data.frame() %>%
left_join(annoInfo) %>%
pull(radiation_field)
## Joining with `by = join_by(sample_id, domain)`
p1 <- Spaniel::spanielPlot_SFE(sfe,
sample_id = "C01",
plot_feat = "domain", ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501264
p2 <- Spaniel::spanielPlot_SFE(sfe,
sample_id = "C01",
plot_feat = "radiation_field", ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501264
p1 + p2
Now that the image annotations are fully loaded in the transcriptomics data can be analysed.
QC metrics about each spot can be calculated using the addPerCellQCMetrics function from scuttle package. By default the total genes per spot are calculated. Additional metrics such as the percentage of mitochondrial genes can be calculated. As this data is from mouse mitochondrial genes start with the prefix “mt-” and can be extracted based on pattern matching.
is_mt <- str_detect(rowData(sfe)$symbol, "^mt-")
sfe <- scuttle::addPerCellQCMetrics(sfe, subsets = list(mito = is_mt))
Spots where no genes are detected can be removed from the remainder of the analysis.
filter <- sfe$detected > 0 & sfe$in_tissue
sfe <- sfe[, filter]
Although we could now progress with the analysis using all samples it is often advisable to analyse each tissue sample or condition separately to begin with. This allows you to identify any section or sample which might be an outlier and should be removed from downstream analysis. For this workshop we example we will split the samples by condition:
sfe_control <- sfe[,sfe$condition == "Control"]
sfe_irradiated <- sfe[,sfe$condition == "irradiated"]
The filtered data can then be then be normalised and dimension reduction steps applied and the spots can be clustered according to transcriptomic similarity. This analysis uses methods taken from the bioconductor Orchestrating Single-Cell Analysis https://bioconductor.org/books/release/OSCA/ work flows. Wrapper function (R/sfe_helper_functions.R) to perform this analysis are implemented in Spaniel for convenience. However for a more comprehensive clustering analysis we recommend following the reading the Scran and Bluster documentation which describes each individual step in the analysis in more detail.
In this example we use a Leidan clustering algorithm as implemented in scran. Leidan is a community detection method. It is important to note that the K parameter is the number of nearest neighbors to use when creating the k nearest neighbor graph. k is related to the resolution of the clustering result, a bigger k will result in lower resolution and vice versa.
sfe_control <- sfe_control %>%
Spaniel::preProcess() %>%
Spaniel::clustCells(K = 100, cluster_function = "leiden")
sfe_irradiated <- sfe_irradiated %>%
Spaniel::preProcess() %>%
Spaniel::clustCells(K = 100, cluster_function = "leiden")
Now we can look at the results of the data both by visualising the clusters on a UMAP plot using plot functions from the scater package and within the spatial context of the tissue section using the spanielPlot function (R/sfe_PlotFunctions.R).
p1_control <- Spaniel::spanielPlot_SFE(sfe_control,
sample_id = "C01",
plot_feat = "domain", ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 500472
p2_control <- Spaniel::spanielPlot_SFE(sfe_control,
sample_id = "C01",
plot_feat = "clust_leiden_K_100",
ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 500472
p1_control + p2_control
scater::plotUMAP(sfe_irradiated, colour_by = "domain")
p1_irradiated <- Spaniel::spanielPlot_SFE(sfe_irradiated,
sample_id = "C06",
plot_feat = "domain", ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501028
p2_irradiated <- Spaniel::spanielPlot_SFE(sfe_irradiated,
sample_id = "C06",
plot_feat = "clust_leiden_K_100",
ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501028
p1_irradiated + p2_irradiated
sfe_control <- sfe_control %>%
Spaniel::clustCells(K = 75, cluster_function = "leiden")
sfe_irradiated <- sfe_irradiated %>%
Spaniel::clustCells(K = 75, cluster_function = "leiden")
p3_control <- Spaniel::spanielPlot_SFE(sfe_control,
sample_id = "C01",
plot_feat = "clust_leiden_K_75",
ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 500472
p3_irradiated <- Spaniel::spanielPlot_SFE(sfe_irradiated,
sample_id = "C06",
plot_feat = "clust_leiden_K_75",
ptSizeMin = 0,
ptSizeMax = 2)
## [1] "got coordinates"
## [1] "metadata"
## [1] "added plot feature"
## [1] "discrete"
## SpatRaster resampled to ncells = 501028
p1_control + p2_control + p3_control
p1_irradiated + p2_irradiated + p3_irradiated
saveRDS(sfe_control, "rObjects/wrapped_sfe_control.rds")
saveRDS(sfe_irradiated, "rObjects/wrapped_sfe_irradiated.rds")